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ABSTRACT 

The  objective  of  this  project  was  to  develop  high-fidelity  models  for  radiation  and  turbulence- 
radiation  interactions  (TRI)  in  high  pressure  combustion  systems  with  extensive  validation  of  the 
aforementioned  models.  New  spectral  radiation  models  for  combustion  gases  at  elevated  pressures, 
for  both  conventional  and  stochastic  solution  methods,  were  developed.  Several  Radiative  Transfer 
Equation  (RTE)  solution  methods  were  newly  designed  (development  of  cell-based  as  well  as 
stochastic  particle-based,  line-by-line  accurate,  photon  Monte  Carlo  methods,  and  formulation 
and  coding  of  higher-order,  3D  spherical  harmonics  as  well  as  simplified  PN  schemes).  To  allow 
simulation  of  high-pressure  laminar  flames  required  modification  of  the  open  source  flow  code 
used  (OpenFOAM)  to  include  differential  diffusion  and  cell-based  stochastic  RTE  solvers.  The 
models  were  validated  by  simulation  of  laminar  high-pressure  hydrogen-air  flames,  as  well  as 
turbulent  methane-air  jet  flames. 

SUMMARY/OVERVIEW 

It  is  well-established  today  that  neglecting  radiation  in  atmospheric  pressure  combustion  systems 
may  lead  to  overprediction  of  temperature  of  up  to  200  °C,  while  using  the  usually-employed 
optically-thin  or  gray  radiation  models  lead  to  underprediction  of  up  to  100 °C  and  more  [1,2]. 
Somewhat  less  well-known  is  the  fact  that,  in  turbulent  combustion  systems,  there  can  be  strong 
interactions  between  turbulence  and  radiation  (TRI).  Probability  density  function  (PDF)  based  cal¬ 
culations  have  shown  that  TRI  always  increases  the  heat  loss  from  a  flame,  and  that  this  additional 
heat  loss  can  reach  60%  of  the  total  and  more,  leading  to  a  reduction  in  the  local  gas  temperature 
of  200  °C  or  more.  This  is  even  more  true  in  high-pressure  combustors  with  their  larger  optical 
thicknesses  [3,4]. 

The  present  research  focused  on  extending  the  Pi’s  previous  work  on  radiation  and  TRI  in  at¬ 
mospheric  pressure  systems  to  the  high  pressures  encountered  in  modern  propulsion  systems.  The 
specific  objectives  were  to  develop,  validate,  and  apply  models  for  RTE  solvers  and  spectral  radi- 


1 


ation  properties  of  the  most  important  combustion  gases  in  high-pressure  combustion,  including 
the  effects  of  TRI.  As  a  result  high-end  models  for  spectral  thermal  radiation  and  TRI  have  been 
advanced  to  a  level  that  is  more  consistent  with  their  importance  in  chemically  reacting  flows. 

TECHNICAL  DISCUSSION 

Four  specific  research  areas  of  this  project  were:  1)  extension  of  spectral  radiation  models  for 
combustion  gases  at  the  elevated  pressures  relevant  to  propulsion  systems,  2)  extension  of  efficient 
and  accurate  RTE  solution  methods  to  high-pressure  flames,  and  3)  perform  a  systematic  and 
detailed  validation  of  the  newly  developed  models. 

Major  Achievements 

(1)  Spectral  radiation  models  for  elevated  pressures  High-fidelity  /.-distribution-based  spectral 
radiation  models  have  been  developed  for  atmospheric  pressure  combustion  [5, 6].  These  models 
have  been  investigated  at  the  elevated  pressures  prevalent  in  propulsion  systems.  In  our  current 
full-spectrum  /^-distribution  (FSK)  model,  k-g  distributions  for  a  mixture  are  constructed  on-the- 
fly  using  a  high-accuracy  database  of  single-species  ^-distributions  together  with  our  narrow-band 
(NB)  mixing  model  [7].  The  NB  mixing  model  [7]  assumes  that  at  the  NB  level  the  absorption 
coefficients  of  gases  are  uncorrelated  (due  to  their  small  overlap)  and,  hence,  the  transmissivities  of 
gases  are  multiplicative.  At  elevated  pressure  much  stronger  line  overlap  is  likely  and  the  accuracy 
of  the  mixing  model  needs  to  be  validated.  We  tested  our  mixing  model  for  pressures  up  to  30  bar 
and  our  findings  are  summarized  in  [8].  Figure  1  shows  the  absolute  errors  {t^xbl  direct- t^lbl  mixing) 
in  NB  transmissivity  calculations  using  our  current  mixing  rule.  Since  important  combustion  gases 
such  as  CO2  and  H20  overlap  primarily  near  2.7  gm.  the  errors  in  transmissivity  calculations  are 
highest  at  this  spectral  location  (as  shown);  everywhere  else  errors  are  negligibly  small.  NB  trans¬ 
missivities  calculated  using  the  mixing  rule  and  single-species  /.'-distribution  databases  (T,hDB  mixing) 
are  essentially  identical  to  LBL-mixing  values  i.e.  T,j  LBL  ,nixing  (not  shown).  Although  our  mixing 
rule  incurs  some  errors  near  3400  cm-1,  due  to  its  extremely  localized  nature  it  has  no  significant 
effect  on  heat  transfer  calculations  (discussed  next). 

A  1-D  homogeneous  isothermal  layer  of  gas  mixture  (1%  C02  and  2%  H20  by  mole)  bounded 
by  cold  black  walls  at  various  pressures  was  considered  keeping  the  pressure  path  length  constant 
at  200  bar.cm.  For  such  a  homogeneous  medium,  the  /.-distribution  method  is  exact  and,  hence, 
is  expected  to  yield  results  within  machine  accuracy  (typically  0.5%).  Any  error  in  addition  to 
machine  error  will  be  due  to  our  mixing  rule  to  construct  mixture  A-z/-cli stri  bution  from  individual 
/.'-^-distributions.  The  nondimensional  heat  fluxes  exiting  the  medium  were  calculated  using  the 
line-by-line  (LBL) ,  the  FSK  and  more  advanced  models,  such  as  the  multi-scale  FSK.  The  results 
are  tabulated  in  Table  1.  It  is  seen  that  errors  in  heat  flux  calculations  using  FSK  and  MSFSK  are 
always  negligibly  small  compared  to  the  LBL  calculations  and  limited  to  1%. 

The  accuracy  of  /.-distribution  methods  was  also  investigated  for  an  inhomogeneous  medium 
at  high  pressure.  A  1-D  medium  containing  a  C02-H20-N2  gas  mixture,  confined  between  cold 
black  walls,  was  considered.  The  mixture  consists  of  two  different  homogeneous  layers  (denoted  as 
left/hot  and  right/cold  layer)  adjacent  to  each  other  with  steps  in  the  temperature  and  concentrations 
of  the  species.  Two  different  total  pressures  of  1  and  30  bar  were  considered.  The  pressure  path 
length  of  the  left/hot  layer  was  kept  constant  at  60  bar-cm  while  the  same  is  varied  for  the  right/cold 
layer.  Such  problems  with  steps  in  species  concentration  and/or  temperature  provide  an  acid  test 
for  these  methods  because  of  their  extreme  inhomogeneity  gradients.  Figure  2  shows  the  results  for 
the  case  where  the  left  layer  contains  10%  CO?,  20%  H20,  and  at  1500  K  whereas  the  right  layer 
contains  20%  CO?,  10%  H20  at  300  K.  It  is  seen  that  the  MSFSK  calculations  are  more  accurate 
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Figure  1:  Comparison  of  NB  transmissivities  calculated  directly  and  using  mixing  rule;  solid  line:  direct, 
dashed  line:  mixing  rule 

(approximately  by  a  factor  of  4)  than  the  single-scale  FSK  results  for  both  pressure.  Accuracy  of 
both  methods  increases  with  increase  in  pressure. 

^-distribution  methods  were  investigated  for  a  more  realistic  but  artificial  methane-air  flame 
(from  Modest  and  Zhang  [9])  at  30  bar  pressure.  Temperature  and  concentration  distributions  for 
C02,  H20  and  CH4  can  be  obtained  from  the  previous  work  of  Modest  and  Zhang  [9].  The  local 
radiative  heat  source  term  is  calculated  using  LBL,  FSK  and  MSFSK  approaches,  employing  the 
P i  method  as  the  RTE  solver,  and  relative  errors  are  determined  by  comparison  with  LBL  as 

N  V  •  <7lbl  -  V  •  <?fsk/msfsk  ,  ,, , 

error(%)  = - - — — - - - - - x  100  (1) 

*  '  <7LBL,max 

The  maximum  error  in  local  radiative  heat  source  term  is  11%  using  LSK  and  2%  using  MSLSK 
(see  Lig.  3)  while  their  atmospheric  counter  parts  have  local  maximum  errors  of  30%  using  LSK 
and  7%  using  MSLSK  (see  Modest  and  Zhang  [9]). 

(2)  Extended  k-distribution  database  for  elevated  pressures  The  original  database  of  Wang  and 
Modest  [10],  which  employed  CDSD-1000  [11]  for  C02  and  HITEMP^S  [12]  for  H20,  was 
updated  to  the  latest  HITEMP-2010  [13],  which  includes  many  million  spectral  lines  activated  at 
high  temperature  and  ensures  accuracy  up  to  4000K.  The  databased  thermodynamic  states  and 
narrow-band  configurations  are  summarized  in  Table  2,  covering  a  wide  range  of  thermodynamic 
states,  larger  than  the  typical  ranges  a  user  may  encounter  in  combustion  simulations.  Lor  exam¬ 
ple,  in  a  low  Mach  number  combustion  simulation,  the  pressure  tends  to  be  almost  constant,  and 
only  a  single  (or  perhaps  very  few)  pressure  values  are  needed.  Similarly,  a  given  combustion 
application  will  encounter  limited  temperature  and  concentration  ranges.  Therefore,  it  is  not  eco¬ 
nomical  to  load  the  entire  database  into  memory,  and  for  this  reason,  the  current  database  utilizes 
dynamic  memory  management,  such  that  each  narrow-band  k-distribution  record  is  only  loaded 
into  memory  at  the  first  time  it  is  requested. 
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Figure  2:  Nondimensional  heat  flux  leaving  an  inhomogeneous  slab  at  total  pressures  of  1  and  30  bar, 
bounded  by  cold  black  walls,  step  changes  in  mole  fraction  and  temperature:  left  layer  contains  10%  CO2, 
20%  H2O,  and  at  1500  K;  right  layer  contains  20%  CO2,  10%  H2O  at  300  K 
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Figure  3:  Relative  error  for  radiative  heat  source  calculations  using  FSK  and  MSFSK  compared  to  LBL  in 
an  artificial  2-D  combustion  chamber  at  30  bar  pressure 
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Table  1:  Nondimensional  Heat  Flux  Exiting  Homogeneous  Isothermal  Medium  at  Various  Pressures  Using 
LBL,  FSK  and  MSFSK  Methods 


T  (K) 

PxF  (bar.cm) 

FBF 

FSK 

MSFSK 

Error  FSK(%) 

Error  MSFSK(%) 

500 

1  bar  x  200  cm 

0.2134 

0.2142 

0.2133 

0.37 

-0.05 

10  bar  x  20  cm 

0.3325 

0.3348 

0.3354 

0.69 

0.87 

20  bar  x  10  cm 

0.3631 

0.3652 

0.3662 

0.58 

0.85 

30  bar  x  6.667  cm 

0.3916 

0.3936 

0.3945 

0.51 

0.76 

1000 

1  bar  x  200  cm 

0.1748 

0.1748 

0.1749 

0.01 

0.07 

10  bar  x  20  cm 

0.2338 

0.2351 

0.2361 

0.56 

0.98 

20  bar  x  10  cm 

0.2474 

0.2483 

0.2499 

0.36 

1.01 

30  bar  x  6.667  cm 

0.2551 

0.2556 

0.2577 

0.20 

1.00 

1500 

1  bar  x  200  cm 

0.1224 

0.1227 

0.1225 

0.27 

0.08 

10  bar  x  20  cm 

0.1486 

0.1485 

0.1491 

-0.07 

0.34 

20  bar  x  10  cm 

0.1538 

0.1537 

0.1551 

-0.70 

0.85 

30  bar  x  6.667  cm 

0.1602 

0.1589 

0.1618 

-0.81 

0.98 

The  accuracy  of  the  narrow-band  database  was  examined  by  repeating  the  sample  calculations 
reported  for  the  previous  version.  Results  are  similar  as  before  and  not  shown  here.  The  relative 
error  is  maintained  to  within  1%  except  for  weak  lines  (mean  absorption  coefficient  below  1  x 
10_6bar_1cm“1),  which  is  consistent  with  previous  observations. 

Differences  between  the  spectroscopy  databases  are  illustrated  in  Fig.  4.  The  pressure -based 
C02  absorption  coefficient  at  a  temperature  of  2000K,  pressure  of  lbar  and  the  zero  mole  fraction 
as  calculated  from  CDSD,  HITEM/J|995  and  HITEMP-2010  are  compared.  The  spectral  range  of 
the  wavenumbers  between  5000  and  5100  cm”1  is  shown,  which  corresponds  to  the  178th  narrow- 
band  of  the  current  database.  The  narrow-band  /.'-distribution  and  the  corresponding  databased 
values  are  also  compared  in  Fig.  5.  For  this  particular  narrow-band,  the  absorption  coefficients 
are  from  HITEMP-2010  are  considerably  larger  than  those  from  CDSD  by  an  approximately  fixed 
value,  and  this  causes  the  narrow-band  /.'-distribution  to  shift  up  by  the  same  amount. 

(3)  Improved  k-distribution  model  for  nonhomogeneous  media  at  elevated  pressures 

Based  on  the  correlation  assumption,  there  are  two  methods  to  find  a  “correlated-k”  value  k* , 

as  shown  in  Fig. 6.  It  is  expected  that  a  reasonable  choice  of  k*  preserves  Planck-mean  absorp¬ 
tion  coefficients  (i.e.,  total  emission),  so  that  radiative  heat  sources  are  accurately  predicted  at  the 
optically-thin  limit.  Figure  7  shows  the  relative  error  in  predicting  pressure -based  Planck-mean 
absorption  coefficients  of  C01  and  H,0  at  the  pressure  of  1  bar  and  zero  concentration  limit  using 
the  same  reference  temperature  of  1000K  but  different  local  temperatures.  For  both  species  con¬ 
sidered,  the  previously  employed  k*  gives  large  errors  as  the  local  temperature  deviates  from  the 
reference  temperature,  while  newly  proposed  k*  gives  relative  errors  mostly  within  1%  even  when 
the  local  temperature  deviates  from  the  reference  temperature  by  more  than  1000K. 

Because  an  error  in  predicting  the  Planck  mean  absorption  coefficient  will  result  in  the  same 
amount  of  error  in  predicting  the  radiative  heat  source  in  optically-thin  conditions  as  encountered 
in  many  laboratory  flames,  the  newly  proposed  k*  is  more  appropriate  and  suggested. 

(4)  Improved  LBL-accurate  photon  Monte  Carlo  method  at  elevated  pressures  In  order  to  use 
a  FBF-accurate  PMC  code,  routines  had  to  be  generated  to  find  statistically  meaningful  emission 
wavelengths  for  photon  bundles  that  are  being  traced.  With  100s  of  millions  of  overlapping  spec¬ 
tral  lines,  finding  efficient  schemes  for  this  purpose  required  challenging  new  models  [14-16].  Re¬ 
cently,  some  improved  wavenumber  selection  schemes  have  been  proposed  for  hypersonic  plasma 
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Table  2:  Precalculated  gas  states  and  narrow-bands 


Parameter 

Sampling 

Number  of  Samples 

Species 

C02,  H20,  CO 

3 

0.1 -0.5:  every  0.1 

5 

0.7: 

1 

Pressure  (bar) 

1.0-14.0:  every  1.0 

14 

15.0-80.0:  every  5 

14 

Total: 

34 

Temperature 

300  -  4000:  every  100 

38 

Mole  Fraction 

0.0  -  1.0:  every  0.25 

5 

200  -  300:  every  10 

10 

300  -  4000:  every  25 

148 

Narrow  Band  (cm-1) 

4000  -  5000:  every  50 

20 

5000-10000:  every  100 

50 

10000-15000:  every  250 

20 

Total: 

248 

3.5rXlCr 


—  CDSD 
HITEMP-2000 

—  HITEMP-2010 


Wavenumber  cm 


Figure  4:  Pressure -based  C02  absorption  coefficients  predicted  by  different  spectroscopy  databases.  Colors 
are  red  (CDSD),  blue  (HITEMP-2000)  and  black  (HITEMP-2010). 
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Figure  5:  Narrow-band  ^-distributions  calculated  from  LBL  absorption  coefficients  and  their  databased  val¬ 
ues.  Lines  are  LBL  calculations,  symbols  are  for  narrow-band  database.  Colors  are  black  (current  version) 
and  red  (previous  version). 


radiation  in  Earth  entry  applications  by  Ozawa  et  al.  [17],  and  by  Feldick  and  Modest  [14].  Unlike 
in  combustion  and  other  high-temperature  systems,  in  Earth  entry  applications,  95%  or  more  of  the 
emitted  energy  comes  from  atomic  species  with  very  few,  but  strong,  electronic  lines.  Contribution 
from  diatomic  species  tend  to  be  minor,  and  diatomic  radiation  is  due  to  relatively  few  rovibrational 
lines.  In  their  scheme,  the  rovibrational  transition  lines  were  separately  databased.  However  there 
are  too  many  rovibrational  transitions  for  combustion  products  such  as  CO2,  H20  and  CO  (several 
hundred  million),  making  it  impractical  to  database  every  individual  line.  By  properly  applying 
the  improved  wavenumber  selection  scheme  to  combustion  gases,  a  hybrid  scheme  was  developed 
for  LBL  photon  Monte  Carlo  simulations  in  order  to  improve  efficiency  as  well  as  to  guarantee 
accuracy. 

In  the  improved  wavenumber  selection  scheme,  emission  from  different  species  is  considered 
separately,  i.e., 


00 

ns  n 

E  f 

7—1  ~ 


Ibrj  9^/- 


(2) 


where  Kr]j  is  the  spectral  absorption  coefficient  of  species  i,  and  1^  is  the  blackbody  intensity. 
Emission  from  individual  species  is  independent  from  another,  i.e.,  there  are  no  overlap  effects 
for  emission.  Rather  than  looking  for  the  appropriate  wavenumber  for  the  mixture,  one  may  first 
determine  the  emitting  species  s  and  then  the  appropriate  wavenumber.  Following  this  idea, 


i=  1  0 

-  +  - 

Etot  Etot 


(3) 


where  Rn  is  a  random  number  and  s  is  the  number  index  of  the  emitting  species.  E:  is  the  total 
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Figure  6:  The  ^-distributions  for  the  evaluation  of  a  “correlated”  k*.  Gas  composition  is  10%  H90  and  90% 
N2,  and  pressure  is  1  bar.  Colors  are  blue  (gas  temperature  of  1000  K)  and  red  (gas  temperature  of  2000K). 
Lines  are  solid  (Planck  function  temperature  1000K)  and  dashed  (Planck  function  temperature  2000K). 


emission  by  species  i  in  the  given  cell,  i.e., 


Ei  =  J  Krhilhll  d/7.  (4) 

0 

In  the  improved  wavenumber  selection  scheme,  first  a  random  number  for  emission  wavenum¬ 
ber,  Rn,  is  drawn,  and  the  emitting  species  5  is  determined  according  to 

j- 1  j 

2  E,  2  Ei 

s  =  j  if  — —  <  R„  <  .  (5) 

J  ns  fl  ns  v  7 

2  Ei  2  Ei 

i=  1  i=  1 

This  ensures  that  the  fraction  Ej/Etot  of  random  numbers  is  employed  to  pick  emission  wavenum¬ 
bers  for  species  i,  in  accordance  with  the  fractional  emission  from  the  species.  Once  the  emitting 
species  is  found,  Rn  is  rescaled  according  to 


RjjEtot  -  2  Ei 

0  <  RrhS  = - - - <  1 

Ej 


(6) 


By  first  determining  the  emitting  species  s  from  Eq.  (5),  an  appropriate  wavenumber  can  be 
found  by  search  through  only  a  single  species  s.  This  improved  scheme,  finding  an  emitting  species 
explicitly  and  and  emission  wavenumber  by  iteration  (in  order  to  eliminate  severe  interpolation 
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T(K) 


T(K) 

Figure  7:  Relative  errors  in  predictions  of  pressure-based  Planck-mean  absorption  coefficients  of  C02  (top) 
and  H-,0  (bottom)  at  different  local  temperatures  using  two  k*  schemes.  is  for  the  original  evaluation 
scheme,  and  k*{  the  newly  developed  method.  The  reference  temperature  is  1000K,  pressure  and  mole 
fraction  are  1  bar  and  zero  concentration  limit  for  both  local  and  reference  states. 
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-^co2=0 . 1 ;  a"h2o=0.1  ;  Aco=0.1 ;  7=650  K 


z  (m) 

Figure  8:  Divergence  of  radiative  flux  for  10%  co2-10%h2o-10%co  mixture  at  650K  using  the  old  scheme 
and  the  hybrid  scheme 

errors)  is  called  “hybrid  scheme.”  In  this  hybrid  scheme  CPU  time  for  the  calculation  of  mixture 
random-number  relations  is  reduced  by  up  to  a  factor  of  10. 

Accuracy  and  efficiency  of  the  old  scheme  and  the  hybrid  scheme  incorporated  into  the  PMC 
model  are  compared  for  a  one-dimensional  example  (in  the  z-direction).  The  example  deals  with 
an  isothermal  gas  at  650K  and  lbar  contained  between  two  parallel  cold  black  walls.  A  mixture 
of  10%C02-10%H20-10%CO  (by  mole  fraction)  and  70%  N2  is  investigated,  and  the  resulting 
divergence  of  the  radiative  flux  is  shown  in  Fig.  8,  comparing  results  from  exact  LBL  calcula¬ 
tions  and  the  old  as  well  as  the  hybrid  emission  wavenumber  selection  schemes.  The  absorption 
coefficient  data  are  taken  from  the  HITEMP2010  database  in  all  three  cases.  It  is  observed  that 
both  Monte  Carlo  schemes  perform  equally  well,  correctly  following  the  LBL  solution.  Using 
10  million  photon  bundles  in  the  PMC  simulation,  the  hybrid  scheme  correctly  follows  the  LBL 
solution,  requiring  3.19  s  to  determine  the  emission  wavenumbers,  as  opposed  to  the  old  scheme, 
which  needed  31.08  s  for  the  identical  example;  i.e.,  the  computational  efficiency  for  wavenumber 
selection  is  improved  about  by  a  factor  of  10. 

(5)  OpenFOAM  Improvements  for  Laminar  Flames  A  laminar-diffusion  combustion  solver, 
taking  into  account  the  effects  of  differential  diffusion,  was  developed  within  the  open  source  CFD 
package  OpenFOAM  [18].  In  addition,  OpenFOAM  was  augmented  to  take  the  effect  of  nonunity 
Lewis  Number  into  account,  which  is  important  in  laminar  flow.  For  the  diffusion  flame  solver 
included  in  OpenFOAM  it  is  assumed  that  the  Lewis  number  is  unity  and  the  diffusivities  of  the 
gas-phase  species  are  equal  to  either  the  mixture  viscosity  or  the  mixture  thermal  diffusivity.  If 
the  assumption  of  unity  Lewis  number  is  relaxed,  an  additional  term  must  be  added  to  the  energy 
equation  due  to  enthalpy  diffusion.  For  slow-flowing  diffusion  flames,  the  absence  of  this  term  can 
result  in  large  errors,  causing  unphysical  temperature  changes  even  when  modeling  the  mixing  of 
two  gases  at  the  same  temperature. 

For  most  chemical  mechanisms,  the  system  of  ODEs  is  stiff  and  ODE  solvers  including  SIBS, 
KRR4  and  RK  in  OpenFOAM  have  been  found  to  be  unstable  for  stiff  systems.  Therefore,  to 
compute  the  chemical  source  terms,  the  CVODE  stiff  ODE  solver  from  the  SunDIALS  package 
[19]  was  incorporated  into  OpenFOAM,  which  is  able  to  solve  both  stiff  and  non-stiff  systems  of 
ODEs. 
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Figure  9:  Two-dimensional  axisymmetric  methane  jet  flame  with  variable  nongray  properties:  (a)  contour 
maps  of  temperature  and  species  mass  fractions;  (b)  radiative  source  at  fixed  axial  location  (x  =  1.0  m)  as 
calculated  by  several  methods. 


(6)  RTE  Solver  Development  To  allow  for  PMC  solutions  in  laminar  flames,  as  well  as  non- 
transported  PDF  turbulent  flame  simulations,  our  stochastic  particle-based  PMC  model  was  mod¬ 
ified  to  allow  for  cell-based  simulations  [15].  The  new  LBL  spectral  module  was  attached  to  both 
the  stochastic-particle  and  the  finite- volume  PMC  codes.  An  example  of  its  calculations  is  given  in 
Fig.  9,  depicting  the  influence  of  radiation  in  a  time-averaged  Sandia  D  type  flow  field,  scaled  up  4 
times  to  make  radiation  important  in  such  an  otherwise  laboratory  scale  flame  (a  snapshot  from  a 
fully  converged  solution  using  our  stochastic  media  LBL  PMC,  and  including  strong  turbulence- 
radiation  interactions).  Figure  9a  shows  temperature  and  concentration  profiles,  indicating  that 
C02  peaks  much  farther  downstream  because  of  considerable  CO  formation  upstream  where  water 
vapor  peaks.  This  puts  great  demands  on  the  spectral  model  (making  absorption  coefficients  “un¬ 
correlated”);  however,  as  seen  from  Fig.  9b,  the  FSK  spectral  model  when  combined  with  our  PMC 
code  gives  essentially  exact  results,  i.e.,  compared  to  LBL-accurate  PMC.  The  strong  gradients  (in 
both  temperature  and  concentrations  in  both  radial  and  axial  directions)  also  challenge  the  simple 
Pi  RTE  solver,  which  incurs  errors  up  to  30%. 

To  improve  the  accuracy  of  conventional  RTE  solvers,  higher-order  PN  methods  were  explored. 
First,  a  simplified  PN,  or  SPN,  model  was  developed,  which  requires  a  set  of  (N  +  l)/2  elliptic  P i 
type  PDEs,  in  the  hope  of  improving  P i  accuracy  at  a  reasonable  price.  The  newly  developed  SPN 
scheme  was  successfully  tested  with  a  number  of  2D  problems  [20],  and  results  are  also  included 
in  Fig.  9b.  It  is  seen  that  maximum  errors  for  SP3  and  SP5  are  reduced  to  about  20%  and  15%, 
respectively,  which  is  slightly  worse  than  would  be  expected  for  the  much  more  complicated  and 
expensive  corresponding  P3  and  P5.  Radiative  heat  sources  for  an  artificial  turbulent  jet  diffusion 
flame  (the  ubiquitous  small  Sandia  Flame  D  scaled  up  4  times  to  make  radiation  more  important), 
employing  different  models,  are  presented  in  Fig.  10.  Results  from  4  different  RTE  solvers  com¬ 
bined  with  3  spectral  models  are  shown.  The  conventional  RTE  solvers  are  paired  with  various 
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k-distribution  spectral  models,  and  P \  also  with  a  LBL  model;  the  PMC  is  always  coupled  with  a 
LBL  spectral  method,  i.e.,  provides  exact  answers  for  comparison.  It  is  observed  that,  when  cou¬ 
pled  with  a  k-distribution,  SP3  eliminates  over  30%  of  the  error  for  P \  compared  PMC,  while  SP5 
recovers  over  one  half.  The  simplified  PN  methods  require  2  or  3  elliptical  equations,  respectively, 
and,  therefore,  provide  relatively  cheap  improvements  over  P\  RTE  evaluations.  However,  when 
coupled  with  the  Planck-mean  gray  model,  P\ ,  SP3  and  SP5  predictions  show  little  difference  to  the 
“optically  thin”  solution  (not  shown),  confirming  that  spectrally  averaged  Planck  mean  absorption 
coefficients  significantly  underestimate  gas  absorption. 

CPU  costs  for  k-distribution  evaluation  of  the  C09  and  H20  mixture  per  cell  is  shown  in  Table 
3.  The  correlation  tables  require  little  cpu  time.  Evaluations  using  the  narrow-band  database 
require  much  more  time  because  of  database  interpolation,  mixing  evaluations  for  each  narrow- 
band,  and  assembly  of  narrow-band  to  full-spectrum  k-distributions.  The  RTE  models,  P 1,  SP3  and 
SP5  require  0.0312s,  0.125s  and  0.172s,  respectively,  for  each  quadrature  on  average.  Because  SP3 
and  SP5  involve  extra  iterations  between  2  and  3  equations  respectively,  their  time  costs  are  more 
than  2  and  3  times  the  P\  cost,  respectively. 


Table  3:  Time  cost  (in  seconds)  for  a  k-distribution  evaluation  per  cell  for  a  C09-H90  mixture. 


Multiplication 

Uncorrelated  mixture 

correlation  tables 

0.0016 

0.006 

narrow-band  database 

0.0549 

0.469 

The  present  calculations  are  based  on  a  mesh  of  3325  cells  and  an  8  point  quadrature  scheme 
for  the  cumulative  k-distribution.  The  total  time  for  k-distribution  assembling  was  5.3  and  1560 
seconds  for  the  correlation  tables  mixed  with  the  multiplication  model  and  for  the  narrow-band 
database  mixed  with  the  uncorrelated  mixture  model,  respectively,  which  represent  computation¬ 
ally  the  cheapest  and  the  most  expensive  k-distribution  models  in  the  present  study.  The  total  CPU 
time  for  the  RTE  solutions  is  0.25,  1  and  1.376  seconds  for  Pi,  SP3  and  SP5,  respectively.  Since 
the  total  computational  time  is  dominated  by  assembling  k-distributions,  simplified  PN  methods 
provide  relatively  cheap  improvements  over  P\  methods. 

Multi-dimensional  PN  implementations  are  mathematically  extremely  complex  [21, 22],  result¬ 
ing  in  a  set  of  N(N  +  l)/2  simultaneous  elliptic  PDEs  with  cross-derivatives.  This  high-order  RTE 
solution  method  has  also  been  implemented  recently,  using  built-in  operator  tools  within  Open- 
FOAM.  These  tools  are  inadequate  for  the  cross-derivatives  required,  as  well  as  for  simultaneous 
solution  of  6,  15,  or  more  PDEs  and,  therefore,  the  present  implementation  is  rather  inefficient.  The 
problem  is  worse  for  the  common  axisymmetric  problems  for  which,  in  theory,  only  ( N  +  l)2/4  (or 
4,  9,  ...)  PDEs  are  needed,  and  which  is  not  exploited  by  the  present  implementation.  Rather  than 
presenting  here  the  extremely  involved  mathematics,  we  show  solutions  for  a  cylindrical  case  to 
demonstrate  PN  solutions  for  axisymmetric  problems.  In  these  types  of  problems,  physical  quan¬ 
tities  such  as  temperature,  heat  flux,  intensity,  chemical  species  concentrations,  etc.,  vary  only 
radially  and  axially,  and  are  considered  as  2-D.  As  a  result,  for  many  of  these  applications,  the 
transport  equations  are  solved  over  a  wedge  in  order  to  reduce  the  computational  effort  in  solv¬ 
ing  these  classes  of  problems.  Here,  a  full  cylinder  rather  than  a  wedge  is  considered  in  order 
to  demonstrate  that  the  resulting  intensity  coefficients  have  azimuthal  dependence  and,  therefore, 
special  care  must  be  taken  regarding  boundary  conditions  for  general  axisymmetric  cases.  For 
this  example  problem,  the  medium  in  the  cylinder  has  strongly  variable  properties  with  an  optical 
thickness  over  a  diameter  at  a  given  height  of  rD  =  4.5.  Results  for  incident  radiation  G  and  nega¬ 
tive  radiative  source  V  •  q  are  shown  in  Figs.  1 1(a)  and  1 1(b).  An  exact  solution  for  this  problem  is 
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Figure  10:  Effects  of  RTE  solvers  on  radiative  heat  flux  divergence  V  •  q  at  three  downstream  locations 
y  =  0.5m  (top),  y  =  1.0m  (middle),  and  y  =  1.4m  (bottom). 
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(a)  Directional  integrated  intensity  (b)  Heat  flux 

Figure  11:  Comparing  incident  radiation  G  and  divergence  of  heat  flux  V  •  q  computed  with  higher  order 
Pn  for  a  cylindrical  enclosure  case.  The  along  the  diameter. 


(a)  /<,  (b)  7°  (c)  if  (d)  l\ 


Figure  12:  Intensity  coefficients  of  P 3  for  cylindrical  problem.  The  resulting  fields  show  that  some  of  the 
intensity  coefficients  exhibit  azimuthal  dependence. 


calculated  directly  from  the  integral  solution  of  the  RTE.  The  solutions  for  P5  and  P1  are  close  to 
the  exact  solution;  however,  increasing  the  order  of  Pn  beyond  N  =  5  changes  the  accuracy  very 
little  for  this  problem.  Surface  plots  of  the  computed  intensity  coefficients  are  shown  in  Fig.  12. 
Values  of  the  intensity  coefficients  pl  and  l\  are  negligible  because  of  the  absent  variations  in  the 
/-direction  and,  therefore,  are  not  shown.  As  this  figure  shows,  the  intensity  coefficients  7?  and  l 2 
have  variations  in  both  the  radial  r  and  polar  6  directions,  whereas  l\\,  has  only  r  dependence  as 
expected.  If  this  problem  were  to  be  applied  on  a  wedge,  new  boundary  conditions  for  7 \  and  722 
must  be  developed. 

While  the  higher  order  P  y - ap p ro x  i  m a t i  o n s  have  been  shown  here  to  be  more  accurate,  they 
typically  require  more  cpu  time.  Details  regarding  the  computational  effort  encountered  on  the  ex- 


Table  4:  CPU  effort  on  a  machine  with  Processor:  Intel  Core  2  Duo  CPU  E8400  3.00GHz  x  2.  The  average 
cpu  time  in  seconds  and  the  number  of  inner  iterations  ( M ). 


Method 

1-D  slab 

60  cells 

2-D  enclosure 
6,889  cells 

Cylinder  Enclosure 
13,520  cells 

Pi 

<0.01 

«0.01 

0.025 

P3 

0.16(8) 

6.22  (24) 

8.87  (10) 

P5 

0.70  (10) 

14.79(16) 

120.48  (44) 

Pi 

2.20(13) 

21.76(10) 

77.45  (12) 

14 


ample  problems  that  were  presented  here  are  displayed  in  Table  4.  This  table  includes  the  number 
of  interior  cells  Nc  for  each  mesh,  the  cpu  time,  and  the  number  of  inner  iterations  (M).  Compared 
with  the  numerical  experiments  from  [23],  where  a  2-D  enclosure  problem  is  solved  with  Nc  =  242 
and  the  cpu  time  was  reported  as  <  Is,  the  application  here  of  P3  on  the  2-D  example  took  6.22s  (an 
order  of  magnitude  greater)  for  a  much  larger  Nc  (also  an  order  of  magnitude  greater).  Comparing 
our  cpu  time  results,  it  appears  that  P3  takes  at  least  2  orders  of  magnitude  longer  than  Pi,  which  is 
considerably  more  than  the  minimum  factor  of  6  (i.e.,  6  simultaneous  PDEs),  due  to  the  presently 
inefficient  OpenFOAM  implementation. 


76.2mm 


25.4mm  .  25.4mm 


2.05mm 

U  1 .525mm 


25.4mm 


Figure  13:  Combustion  chamber  geometry. 

(7)  High-Pressure  Flame  simulations  We  studied  a  few  high-pressure,  laminar  hydrogen-air 
flames  numerically,  similar  to  flames,  which  can  be  studied  experimentally  at  the  Edwards  AFB 
EC-1  facility  in  future.  These  flames  were  studied  using  both  conventional  RTE  solvers  (Pi, 
SP3  and  SP5,  together  with  the  FSK  spectral  model),  as  well  as  LBL-accurate  PMC  calcula¬ 
tions  [15, 24].  The  modeled  chamber  is  an  axisymmetric  version  of  the  EC-1  facility  and  is  shown 
in  Fig.  13.  The  chamber  has  a  radius  of  25.4mm  with  an  exit  radius  of  12.7mm.  The  central  jet  has 
a  diameter  of  3.05mm,  and  an  annular  jet  surrounding  it  has  an  inner  diameter  of  4.1mm.  The  wall 
thickness  of  the  inner  jet  is  neglected.  The  central  jet  supplies  air  at  a  temperature  of  300K  and 
a  mass  flow  rate  of  0.143g/s.  The  annular  jet  supplies  hydrogen  at  a  temperature  of  300K  and  a 
mass  flow  rate  of  4.16  x  10~3g/s.  Simulations  were  performed  for  three  pressures,  viz.,  1,  5  and  30 
bar.  In  all  three  cases,  the  same  mass  flow  rates  are  used.  Because  gas  dynamic  viscosity  has  little 
dependency  on  pressure,  the  Reynolds  number  based  on  air  mass  flow  rate,  viscosity  and  inner  jet 
diameter  is 

4m 

Re  = - =  3268  (7) 

nD/U 

for  all  pressures.  This  Reynolds  number  is  much  lower  than  the  transient  Reynolds  number  of 
laboratory  jet  flames,  such  as  Sandia  Flames  [25]  and  is  expected  to  result  in  a  laminar  flame. 
The  wall  temperature  is  set  to  400K.  A  two-dimensional  cylindrical  grid  with  640  cells  is  used  to 
discretize  all  model  equations,  and  is  refined  near  the  flame  front.  Sample  results  are  shown  in 
Fig.  14,  showing  temperature  maps  for  the  3  different  pressures,  as  calculated  without  radiation 
and  with  radiation,  using  the  simplest  Pi/FSK  solver.  Without  radiation,  maximum  temperature 
increases  with  pressure,  but  at  the  same  time  the  high-temperature  region  becomes  narrower  be¬ 
cause  of  faster  reaction  rates.  Radiative  cooling  is  evident  even  for  the  optically  thin  lbar  flame, 
becoming  stronger  at  5bar,  and  overwhelming  at  30bar.  This  is  better  seen  by  looking  at  a  single 
cross-section  (taken  to  be  the  chamber  exit,  giving  an  indication  of  the  total  heat  loss  from  the 
flame),  as  shown  in  Fig.  15:  radiation  cools  down  the  lbar  flame  by  roughly  100K  (except  near 
the  centerline  where  one  sees  a  temperature  increase,  caused  by  downstream  absorption),  the  5bar 
flame  by  200K,  and  the  30bar  flame  by  600K!  At  lbar  the  flame  is  optically  thin,  and  all  RTE 
solvers  return  the  same  result  (since  there  is  little  self-absorption).  Even  at  higher  pressures  the 
differences  in  results  between  Py  and  higher-order  solvers  are  small,  indicating  that  P\  returns 
accurate  results. 
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Figure  14:  Temperature  profiles  for  laminar  model  combustor  at  3  pressures;  left  frames:  no  radiation,  right 
frames:  radiation  calculated  by  P  \  /FS K  model 
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Figure  15:  Temperature  profiles  at  the  exit  cross-section,  x  =  0.1 2m 
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